##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(RColorBrewer)
library(ggpubr)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--ccf_file"), type = "character") ,
    make_option(c("--gene_list"), type = "character") ,
    make_option(c("--ith_file"), type = "character") ,
    make_option(c("--ith_sample_file"), type = "character") ,
    make_option(c("--sample_info"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    ccf_file <- "~/20220915_gastric_multiple/dna_combinePublic/mutationTime/result/All_CCF_mutTime.addShare.tsv"
    sample_info <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.list"
    gene_list <- "~/20220915_gastric_multiple/dna_combinePublic/images/selectGCClone/GCClone_gene.all_record.list"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/images/ITH"
    ith_file <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/Evolution_Mode/ITH.compute.uniqueNormal.tsv"
    ith_sample_file <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/Evolution_Mode/ITH.compute.allSample.tsv"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gene_list <- opt$gene_list
ccf_file <- opt$ccf_file
sample_info <- opt$sample_info
out_path <- opt$out_path
ith_file <- opt$ith_file
ith_sample_file <- opt$ith_sample_file

###########################################################################################

dir.create(out_path , recursive = T)
col <- c( "#006699","#DDA520"  )

###########################################################################################

dat_ccf <- fread( ccf_file )

dat_info <- data.frame(fread(sample_info))
dat_gene <- data.frame(fread(gene_list , header = T))
dat_ith <- data.frame(fread(ith_file))
dat_ith_sample <- data.frame(fread(ith_sample_file))

colnames(dat_gene) <- "Gene_Symbol"

###########################################################################################
## 排除癌前病变的驱动基因
## 当癌前和癌共享这个基因时，其ITH仍然很高，说明非关键的驱动因素
## dat_gene <- subset(dat_gene , !(Gene_Symbol %in% dat_pre$Gene_Symbol))

###########################################################################################

col<-brewer.pal(12,"Set3")
col_point <- brewer.pal(9,"Set1")
Variant_Type <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

dat_ccf$Location <- paste( dat_ccf$Chr , dat_ccf$Start_Position , 
    dat_ccf$REF , dat_ccf$ALT , sep=":"  )

###########################################################################################
judgeDriver <- function(dat_info = dat_info , dat_ccf = dat_ccf){
    result_driver <- c()

    for( Sample in unique(dat_info$ID) ){

        tumors <- subset( dat_info , ID == Sample )$Tumor

        ## 确定driver突变
        tmp <- subset( dat_ccf , Sample %in% tumors & Variant_Classification %in% Variant_Type & Hugo_Symbol %in% dat_gene$Gene_Symbol )
        #tmp <- subset( tmp , Class == "IM" | (Class!="IM" & CCF_adj >= 0.6) )

        ## 判断突变的数量
        tmp_driver <- tmp %>% 
        group_by(Location) %>%
        summarize( MutTumor = length(Sample) )

        if( nrow(tmp_driver) == 0){
            class <- "NoDriver"
        }else if( nrow(tmp_driver) > 0 ){
            share_driver <- length(which(tmp_driver$MutTumor == length(tumors)))

            ## 判断是否存在部分IM与GC共享
            tmp1 <- tmp
            tmp_driver1 <- tmp1 %>% 
            group_by(Location ) %>%
            summarize( MutClass = paste0(Class , collapse = "_") )
            share_driver1 <- length(which(tmp_driver1$MutClass %in% c("IM_IGC" , "IM_DGC" , "IGC_IM" , "DGC_IM")))

            if(share_driver > 0 ){
                class <- "ShareDriver"
            }else if(share_driver1 > 0){
                class <- "PartDriver"
            }else{
                class <- "PrivateDriver"
            }
        }

        tmp_res <- data.frame( Normal = Sample , DriverClass = class )

        result_driver <- rbind( result_driver , tmp_res )
    }

    return(result_driver)
}

###########################################################################################
## 判断肿瘤的类型
## 1、无驱动突变
## 2、已报道的driver突变出现在Trunk
## 3、driver突变均在分支
dat_info_IGC_DGC <- subset( dat_info , Type == "IM + IGC + DGC" )
dat_info <- subset( dat_info , Type != "IM + IGC + DGC" )

result_driver <- judgeDriver(dat_info = dat_info , dat_ccf = dat_ccf)
result <- merge( dat_ith , result_driver , by.x = "ID" , by.y = "Normal" )

result_driver_IGC_DGC_IGC <- judgeDriver(dat_info = subset( dat_info_IGC_DGC , Class != "DGC" ) , dat_ccf = dat_ccf)
result_driver_IGC_DGC_DGC <- judgeDriver(dat_info = subset( dat_info_IGC_DGC , Class != "IGC" ) , dat_ccf = dat_ccf)
result_IGC <- merge( subset( dat_ith , Type=="IM + IGC + DGC(IGC)") , result_driver_IGC_DGC_IGC , by.x = "ID" , by.y = "Normal" )
result_DGC <- merge( subset( dat_ith , Type=="IM + IGC + DGC(DGC)") , result_driver_IGC_DGC_DGC , by.x = "ID" , by.y = "Normal" )

###########################################################################################
## 合并异质性和驱动突变
result <- rbind( result , result_IGC , result_DGC )

result$DriverClass_combine <- ifelse( result$DriverClass=="ShareDriver" , "TrunkDriver" , "NoTrunkDriver")
result$DriverClass_combine <- factor( result$DriverClass_combine , levels = c( "TrunkDriver" , "NoTrunkDriver") , order = T )
result_final <- result

###########################################################################################

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("p == ",down," %*% 10","^",up)
    return(text)
}

###########################################################################################
## 计算两者的p值
p <- wilcox.test(subset(result , DriverClass_combine=="TrunkDriver")$ITH , subset(result , DriverClass_combine=="NoTrunkDriver")$ITH)$p.value

if( p < 0.01 ){
    p_text <- trans(p)
}else{
    p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
}

result$p_text <- ""
result$p_text[1] <- p_text

if(1!=1){
    out_name <- paste0(out_path , "/ITH.combineNormal.TrunkDriver.pdf")
    plot <- ggplot(data=result,mapping = aes(x=DriverClass_combine,y=ITH))+
      geom_boxplot(lwd=1.5,aes(color=DriverClass_combine) , outlier.colour = NA) +
      geom_jitter(position=position_jitter(0.2),aes(color=DriverClass_combine)) +
      scale_color_npg() +
      #facet_grid(.~Type)+
      xlab(NULL) +
      ylab('ITH')+
      theme_bw() +
      ylim(0.2,1) +
      geom_text(aes(label=p_text , y = 1 ,x = 1.5),parse = TRUE,size=4)+
      theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='none',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            legend.text = element_text(size = 8,color="black",face='bold'),
            axis.text.x = element_text(size = 10,color="black",face='bold'),
            axis.text.y = element_text(size = 10,color="black",face='bold'),
            axis.title.x = element_text(size = 10,color="black",face='bold'),
            axis.title.y = element_text(size = 12,color="black",face='bold'),
            axis.line = element_line(size = 0.5)) 
    ggsave(file=out_name,plot=plot,width=3,height=5)
}

###########################################################################################
## 在IGC和DGC两者中看TurnkDriver
result$p_sub_text <- ""
result$Type <- ifelse( result$Type == "IM + IGC + DGC(IGC)" , "IM + IGC" , result$Type )
result$Type <- ifelse( result$Type == "IM + IGC + DGC(DGC)" , "IM + DGC" , result$Type )

p <- wilcox.test(
    subset(result , Type == "IM + IGC" & DriverClass_combine=="TrunkDriver")$ITH , 
    subset(result , Type == "IM + IGC" & DriverClass_combine=="NoTrunkDriver")$ITH)$p.value
if( p < 0.01 ){
    p_text <- trans(p)
}else{
    p_text <- paste0( "p == " , round(as.numeric(p) , 4) ) 
}
result[which(result$Type == "IM + IGC"),]$p_sub_text[1] <- p_text

p <- wilcox.test(
    subset(result , Type == "IM + DGC" & DriverClass_combine=="TrunkDriver")$ITH , 
    subset(result , Type == "IM + DGC" & DriverClass_combine=="NoTrunkDriver")$ITH)$p.value
if( p < 0.01 ){
    p_text <- trans(p)
}else{
    p_text <- paste0( "p == " , round(as.numeric(p) , 4) ) 
}
result[which(result$Type == "IM + DGC"),]$p_sub_text[1] <- p_text

result$Type <- factor(result$Type , levels = c("IM + IGC" , "IM + DGC") )

col_tmp <- c(
        rgb(red=247,green=184,blue=71,alpha=255,max=255) ,
        rgb(red=2,green=100,blue=190,alpha=255,max=255) 
    )

out_name <- paste0(out_path , "/ITH.combineNormal.IGC_DGC.TrunkDriver.pdf")
plot <- ggplot(data=result,mapping = aes(x=Type,y=ITH,color = DriverClass_combine ))+
    geom_boxplot(alpha =1 , size = 0.9 , width = 0.6 , outlier.shape = NA) +
    geom_point(position = position_jitterdodge() , aes(color=DriverClass_combine)) +
    scale_color_manual(values=col_tmp) +
    xlab(NULL) +
    ylab('Heterogeneity index')+
    theme_bw() +
    ylim(0.2,1.1) +
    geom_text(aes(label=p_text , y = 1.1 ,x = 1.5),parse = TRUE,size=4,color = "black")+
    geom_text(aes(label=p_sub_text , y = 1 ,x = Type),parse = TRUE,size=3,color = "black")+
    theme(
        legend.position = 'bottom',
        legend.title = element_blank() ,
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        panel.background = element_blank(),
        panel.border = element_blank(),
        plot.title = element_text(size = 12,color="black",face='bold'),
        legend.text = element_text(size = 8,color="black",face='bold'),
        axis.text.y = element_text(size = 12,color="black",face='bold'),
        axis.title.x = element_text(size = 12,color="black",face='bold'),
        axis.title.y = element_text(size = 12,color="black",face='bold'),
        axis.text.x = element_text(size = 12,color="black",face='bold') ,
        axis.ticks.length = unit(0.2, "cm") ,
        strip.text.x = element_text(size = 15, colour = "black",face='bold') ,
        axis.line = element_line(size = 0.5)
    ) 

ggsave(file=out_name,plot=plot,width=3,height=4)

out_name <- paste0(out_path , "/ITH.combineNormal.IGC_DGC.TrunkDriver.tsv")
write.table( result , out_name , row.names = F , sep = "\t" , quote = F )

## 判断每个人是否共享驱动突变
result_normal <- data.frame(Normal = result$ID , DriverClass = result$DriverClass , DriverClass_combine = result$DriverClass_combine)
result_normal <- unique(result_normal)
out_name <- paste0(out_path , "/DriverClass.tsv")
write.table( result_normal , out_name , row.names = F , sep = "\t" , quote = F )

